This R code includes the analyses and figures of the SMART manuscript. This R code includes analyses of Nora Liebers, Peter-Martin Bruch, Miguel Hernandez, Junyan Lu and Tobias Terzer.

1 1 Setup

library(knitr)
options(digits=3, width=80)
opts_chunk$set(echo=TRUE,tidy=FALSE,include=TRUE)

1.1 1.1 Loading libraries

library(ggpmisc)
suppressMessages(library(RColorBrewer))
suppressMessages(library(ggrepel))
suppressMessages(library(gghighlight))
suppressMessages(library(ggbeeswarm))
suppressMessages(library(cowplot))
suppressMessages(library(gridExtra))
suppressMessages(library(drc))
suppressMessages(library(tidyverse))
detach("package:dplyr") 
suppressMessages(library(dplyr))
library("survminer")
library("survival")
library("maxstat")
library("forestmodel")
library("forestplot")
library("patchwork")
library("table1")
library("glmnet")
library("corrplot")
library("writexl")

#maintain function of specific R packages
rename <- dplyr::rename
count  <- dplyr::count
filter <- dplyr::filter

1.2 1.2 Define figure themes

fontsize=18
fontsize_small=10
starsize <- 5
statsize <- 5


theme_figures<- theme_classic()+
  theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
        axis.text.x =element_text(size=fontsize),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        axis.line = element_line(size=0.6),
        axis.ticks = element_line(size=1.5),
        legend.text = element_text(size=fontsize),
        legend.title = element_text(size=fontsize+2))

theme_figures_wo_legend<- theme_classic()+
  theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
        axis.text.x =element_text(size=fontsize),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        axis.line = element_line(size=0.6),
        axis.ticks = element_line(size=1.5),
        legend.position = "none")

theme_figures_x_angle_45<-theme_classic()+
  theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
         axis.text.x =element_text(size=fontsize, angle=45, vjust=1, hjust=1),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        axis.line = element_line(size=0.6),
        axis.ticks = element_line(size=1.5),
        legend.text = element_text(size=fontsize),
        legend.title = element_text(size=fontsize+2))
  
theme_figures_x_angle_45_wo_legend<-theme_classic()+
  theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
        axis.text.x =element_text(size=fontsize, angle=45, vjust=1, hjust=1),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        axis.line = element_line(size=0.6),
        axis.ticks = element_line(size=1.5),
        legend.position = "none")
  
 theme_figures_x_angle<-theme_classic()+
  theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
         axis.text.x =element_text(size=fontsize, angle=90, vjust=0.5, hjust=1),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        axis.line = element_line(size=0.6),
        axis.ticks = element_line(size=1.5),
        legend.text = element_text(size=fontsize),
        legend.title = element_text(size=fontsize+2))

theme_figures_x_angle_wo_legend<-theme_classic()+
  theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
         axis.text.x =element_text(size=fontsize, angle=90, vjust=0.5, hjust=1),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        axis.line = element_line(size=0.6),
        axis.ticks = element_line(size=1.5),
        legend.position = "none")

# used colors in combination with alpha=0.85
## primary colors
red1<-"orangered3"
gray1<-"slategray4"
blue1<-"deepskyblue4"
 
## secondary colors
red2<-"orangered4"
blue2<-"deepskyblue3"

1.3 1.3 Define functions

1.3.1 1.3.1 Functions for Table 1

my.render.cont <- function(x) {
    with(stats.apply.rounding(stats.default(x), digits=3), c("",
        "Median [Min, Max]"=sprintf("%s [%s,  %s]", MEDIAN, MIN, MAX)))
}

my.render.cont_IQR <- function(x) {
    with(stats.apply.rounding(stats.default(x), digits=3), c("",
        "Median [Min, Max]"=sprintf("%s [%s,  %s]", MEDIAN, IQR)))
}

my.render.cat <- function(x) {
    c("", sapply(stats.default(x), function(y) with( 
      y,sprintf("%d (%0.0f %%)", FREQ, PCT))))
}

1.3.2 1.3.2 Functions for Analyses

#standard or Welch's two-sample t-test 
t_test <- function(yourTab, var_eq=T){
  pVal <- t.test(metric ~ response, yourTab, var.equal = var_eq)$p.value
}

m_reg <- function(yourTab){
  #pvalue from the response comparison after accounting for tumor infiltration
  return( summary( lm(AUC ~ response + t_infiltration, data = yourTab) )$coefficients[2,4] )
}


###Dose-response curves
fit_log_curve <- function(df, drugname){
  #Fits a 5-parameter logistic model for each in-vivo response and returns:
  # 1) coefficients (b,c,d,e,f) for each model.
  # 2) r^2 between predicted and real values for each model.
  # 3) a df with 1k predicted viability values 
  #    for each in-vivo response. This can then be used for building a smooth dose-response curve. 
  
  fun1 <- function(df, drugname){
    # Fit individual curves (5 parameter logistic worked best)
    # parameters and meanings
    # b, c, d, e, f
    # slope, lower, upper, ED50, asymmetry
    mR <- drm(normVal ~ concentration, data = df[df$response == "R",], fct = LL2.5(), 
                type = "continuous", upperl = c(NA, 1, NA, 1, NA), lowerl = c(NA, 0, NA, 0, NA)) 
    mPD <- drm(normVal ~ concentration, data = df[df$response == "PD",], fct = LL2.5(),
                 type = "continuous", upperl = c(NA, 1, NA, 1, NA), lowerl = c(NA, 0, NA, 0, NA)) 
    #get param coefficients
    coefR <- mR$coef 
    coefPD <- mPD$coef 
    #evaluate fit
    r2_mR <- round( cor(predict(mR),df[df$response == "R",]$normVal, use = "pairwise.complete.obs")^2, 2)
    r2_mPD <- round( cor(predict(mPD),df[df$response == "PD",]$normVal, use = "pairwise.complete.obs")^2, 2)
    #predict 1k concentrations to smoothen curves in line plot (dose-response curve)
    df_conc <- data.frame(concentration = seq(min(df$concentration), max(df$concentration), length.out=1000))
    predi_R <- predict(mR, newdata = df_conc )
    predi_PD <- predict(mPD, newdata = df_conc )
      
    return( list( coefR, 
                  coefPD, 
                  r2_mR, 
                  r2_mPD,
                  tibble(name = rep(drugname, 2000), 
                         response = c(rep("R", 1000), rep("PD", 1000)) , 
                         concentration = c( seq(min(df$concentration), max(df$concentration), length.out=1000),
                                            seq(min(df$concentration), max(df$concentration), length.out=1000) ) , 
                         viab_pred = c(predi_R,predi_PD)) ) )
  }
  #function if error
  errorfun <- function(){
      #if fitting error, return a df with NAs
      return( list( c(NA,NA,NA,NA,NA), 
                    c(NA,NA,NA,NA,NA),
                    NA,
                    NA,
                    tibble(name = drugname, response = NA, concentration = NA, viab_pred = NA) ) )}
  
  #call functions
  out_lst <- tryCatch( fun1(df, drugname), 
                       error = function(e){T} )
  #handle error and return
  if ( typeof(out_lst) == typeof(T) ){return( errorfun() )} else{return( out_lst )}
}

#Junyan's 
#To calculate the area under the curve based on the trapezoidal rule
calcAUC <- function(value, conc) {
    #value is a vector of normalized viability values
    #conc is a vector of concentrations.
    
    valueConc <- data.frame(viab=value, conc=conc)
    valueConc <- valueConc[order(valueConc$conc),]
    areaTotal <- 0
    for (i in seq(length(value)-1)) {
        areaEach <- (valueConc$viab[i] + valueConc$viab[i+1])*log(valueConc$conc[i+1]/valueConc$conc[i])*0.5
        areaTotal <- areaTotal + areaEach
    }
    areaTotal/log(valueConc[nrow(valueConc),]$conc/valueConc[1,]$conc)
}

#Junyan's miscellaneous custom functions: library(jyluMisc)
#it creates a pdf doc where it arranges the plots contained in a plotList on a grid 
makepdf <- function(x, name, ncol = 3, nrow = 2, figNum =NULL, width =20, height = 12) {
  require(gridExtra)
  #figs per page
  if (is.null(figNum)) figNum = ncol*nrow
  if (length(x) == 0) return(NULL)
  #create empty pdf doc with dims
  pdf(name, width = width, height = height)
  #iterate over len of plots in plotList by figNum (e.g., 8 by 8 = 1, 9, 17...)
  for (i in seq(1, length(x), by = figNum)) {
    #print(names(x)[i])
    j  <- min(i + figNum - 1, length(x))
    #arrange plots in grid 
    do.call(grid.arrange, c(x[i:j], ncol = ncol, nrow = nrow))
  }
  dev.off() %>% invisible #close pdf
}


# Functions for elastic net regression #
repeated_glmnet <- function(X, Y, alpha = alpha, type.measure = "auc", nlambda, nfolds, ncv = 1000, family = "binomial"){
  folds.list <- lapply(1:ncv, function(x){return(balanced_folds(Y, nfolds))})
  lasso.models <- lapply(X = folds.list, FUN = function(folds){
    return(cv.glmnet(X, Y, alpha = alpha, type.measure = type.measure, intercept = TRUE, standardize = FALSE, nlambda = nlambda, foldid = folds, family = family))
  })
  
  tmp <- lapply(lasso.models, FUN = function(x){
    lambda.min <- x$lambda.min
    index.min <- x$lambda == lambda.min
    cvm <- x$cvm[index.min]
    coef <- x$glmnet.fit$beta[,index.min]
    names(coef) <- rownames(x$glmnet.fit$beta)
    return(list(lambda = lambda.min, cvm = cvm, coef = coef))
  })
  
  coefficients <- matrix(unlist(lapply(tmp, FUN = function(x){
    return(x$coef)})), ncol = length(tmp), byrow = FALSE)
  rownames(coefficients) <- names(tmp[[1]]$coef)
  
  cvm <- sapply(tmp, FUN = function(x){
    return(x$cvm)
  })
  
  lambda <- sapply(tmp, FUN = function(x){
    return(x$lambda)
  })
  return(list(coef = coefficients, cvm = cvm, lambda = lambda))
}

give_repeated_lasso_results <- function(lasso.models, type.measure = "auc"){
  coef_mat <-lasso.models$coef
  median_coef <- apply(coef_mat, MARGIN = 1, FUN = function(x){
    return(ifelse(any(x != 0), median(x[x != 0]), 0))})
  
  proportion_selected <- apply(coef_mat, MARGIN = 1, FUN = function(x){
    return(sum(x != 0)/length(x))})
  cvms <- lasso.models$cvm
  
  median_auc <- sapply(1:nrow(coef_mat), FUN = function(cur_index){
    relevant_cols <- which(coef_mat[cur_index,] != 0)
    return(ifelse(length(relevant_cols) > 0, median(cvms[relevant_cols]), 0.5))})
  
  tmp.frame <- data.frame(Covariate_label = rownames(coef_mat), medianOR = median_coef, Proportion_selected = proportion_selected, median_AUROC  = median_auc)
  rownames(tmp.frame) <- tmp.frame$Covariate_label
  tmp.frame$medianOR <- case_when(tmp.frame$Covariate_label %in% c("Vidaza", "Azaguanin", chemotherapies) ~ exp(tmp.frame$medianOR/10), TRUE ~ exp(tmp.frame$medianOR))
  tmp.frame <- tmp.frame[order(tmp.frame$Proportion_selected, decreasing = TRUE),]
  rownames(tmp.frame)
  
  #Compute proportions matrix
  coef_mat <- coef_mat[rownames(tmp.frame),]
  mat.tmp <-matrix(0, nrow = nrow(coef_mat), ncol = nrow(coef_mat))
  for(i in 1:nrow(mat.tmp))
  {
    relevant_cols <- which(coef_mat[i,] != 0)
    for(j in 1:nrow(mat.tmp))
    {
      if(length(relevant_cols) != 0){
        mat.tmp[i,j] <- sum(coef_mat[j, relevant_cols] != 0) /length(relevant_cols)
      }
      else
      {
        mat.tmp[i,j] <- 0
      }
    }
  }
  rownames(mat.tmp) <- colnames(mat.tmp) <- rownames(coef_mat)
  diag(mat.tmp) <- proportion_selected[rownames(mat.tmp)]
  return(list(overview_frame = tmp.frame, proportions.mat = mat.tmp))
}

# function to select folds in cv in a balanced fashion
balanced.folds.new <- function (y, nfolds = min(min(table(y)), 10)) 
{
  totals <- table(y)
  fmax <- max(totals)
  nfolds <- min(nfolds, fmax)
  nfolds = max(nfolds, 2)
  folds <- as.list(seq(nfolds))
  yids <- split(seq(y), y)
  bigmat <- matrix(NA, ceiling(fmax/nfolds) * nfolds, length(totals))
  for (i in seq(totals)) {
    # cat(i)
    if (length(yids[[i]]) > 1) {
      bigmat[seq(totals[i]), i] <- sample(yids[[i]])
    }
    if (length(yids[[i]]) == 1) {
      bigmat[seq(totals[i]), i] <- yids[[i]]
    }
  }
  smallmat <- matrix(bigmat, nrow = nfolds)
  smallmat <- pamr:::permute.rows(t(smallmat))
  res <- vector("list", nfolds)
  for (j in 1:nfolds) {
    jj <- !is.na(smallmat[, j])
    res[[j]] <- smallmat[jj, j]
  }
  return(res)
}

balanced_folds <- function (class.column.factor, cross.outer){
  sampleOfFolds <- get("balanced.folds.new")(class.column.factor, nfolds = cross.outer)
  permutated.cut <- rep(0, length(class.column.factor))
  for (sample in 1:cross.outer){
    # cat(sample, "\n")
    permutated.cut[sampleOfFolds[[sample]]] <- sample
  }
  return(permutated.cut)
}

2 2 Data read in and preprocessing

2.1 2.1 File path

# This filepath can be modified if necessary

filepath<- "../../"

2.2 2.2 Load processed data

load(paste0(filepath,"data/SMARTrial_data.RData"))

screenData<-left_join(screenData, df_drugs, by="Drug")
data_all<-left_join(data_all, df_drugs, by="Drug")

Explanation of all datasets:
- screenData: processed ex-vivo raw data
- data_all: includes in-vivo and ex-vivo data of all 80 patients after data quality control
- df_chemo: includes only patients with chemotherapy and evaluable response (n=46) after data quality control, includes AUC and normalized viability per concentration and drug
- df_drugs: list and metadata of drugs used ex-vivo

3 3 Overview drugs

3.1 3.1 Figure: overview drugs

# prepare data
drugs_overview <-df_drugs %>% dplyr::count(pathway, FDA_approved) %>% 
         left_join(., dplyr::count(df_drugs, pathway), by="pathway")

drugs_overview$pathway <-factor(drugs_overview$pathway, levels=c("ABL (BCR)", "ALK", "Apoptosis", "BCR", "Bromodomain", "CDK", "Cell cycle" , "Chemotherapy", "DDR", "FLT3", "Hedgehog", "HGF", "Histone acetyltransferase", "Histone deacetylase", "Histone methyltransferase", "IDH", "JAK/STAT", "MAPK", "Notch", "Nuclear export", "PI3K/AKT/mTOR", "PKC", "PLK", "Promiscuous" , "Proteasome"    , "Stress response", "TLR", "TNF/NFKB" ), 
labels=c("BCR/ABL", "ALK", "Apoptosis", "B cell receptor", "Bromodomain", "Cyclin dependent kinases", "Cell cycle control" , "Chemotherapy", "DNA damage response", "FLT3", "Hedgehog", "HGF", "Histone acetyltransferase", "Histone deacetylase", "Histone methyltransferase", "IDH", "JAK/STAT", "MAPK", "Notch", "Nuclear export", "PI3K/AKT/mTOR", "PKC", "PLK", "Other" , "Proteasome", "Stress response", "TLR", "TNF/NFKB" ))

# define figure theme
theme_tmp<-theme_classic()+
  theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
         axis.text.x =element_text(size=fontsize, angle=45, vjust=1, hjust=1),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        axis.line = element_line(size=0.6),
        axis.ticks = element_line(size=1.5),
        legend.text = element_text(size=fontsize),
        legend.title = element_text(size=fontsize+2),
        legend.position = c(0.7, 0.9))
  
# plot
ggplot(data=drugs_overview, aes(x=reorder(pathway, -n.y), y=n.x, fill=as.character(FDA_approved)))+
         geom_bar(stat = "identity", width=0.6)+
         scale_fill_manual(values=c("#cbd1d5", "#f2a4a3"), name="Drug type", label=c("Clinical development/tool compound", "FDA approved"))+
         ylab("\n\nNumber of drugs\n")+
         xlab("")+
         theme_tmp

# ggsave(width = 12, height =8, path=(paste0(filepath,"04_Figures")), filename = "Figure_1B_Compound_library.svg")

remove(drugs_overview, theme_tmp)

3.2 3.2 Supplementary table: Drug list

For appendix: includes the pathway, target, Supplier of each drug

df_drugs_supp<-df_drugs%>%select(Drug, pathway, target, Supplier)%>%filter(Drug!="Bortezomib + Cytarabine(800nM)")%>%rename(Compound="Drug")%>%rename(Main.targets_mode.of.action="target")%>%rename(Simplified.drug.class="pathway")

df_drugs_supp %>% DT::datatable(width = 700)
# write_xlsx(df_drugs_supp,"04_Figures/table_S1_drugs.xlsx")

remove(df_drugs_supp)

4 4 Patient characteristics

4.1 4.1 Table 1

# preparation
df_invivo$diagnosis_summarized <-factor(df_invivo$diagnosis, levels=c("AML", "ALL","DLBCL", "BL", "FL", "MCL", "CLL", "B-PLL", "T-PLL", "T-NHL" ), 
                       labels=c("AML", "ALL", "DLBCL", "Burkitt lymphoma", "Follicular lymphoma", "Mantle cell lymphoma", "CLL", "B-PLL", "T-PLL", "T-cell lymphoma, NOS" ))

label(df_invivo$age) <-"Age (years)"
label(df_invivo$diagnosis_summarized) <-"Diagnosis"
label(df_invivo$time_from_ED_dich) <-"Time from initial diagnosis (months)"
label(df_invivo$lines_pretreatment_cat) <-"Prior lines of treatment"
label(df_invivo$tumor_infiltration) <-"Tumor infiltration of sample"
label(df_invivo$Treatment_type) <-"Prescribed treatment at study entry"

df_invivo$material <-factor(df_invivo$material, levels=c("PB", "BM", "LN"),
                                       labels=c("Periphal blood", "Bone marrow", "Lymph node"))
label(df_invivo$material) <-"Tumor sample origin"

df_invivo<-df_invivo %>% 
  mutate(Treatment_type=gsub("mall molecules","mall molecule inhibitors",Treatment_type))

# table 1 printing
table1(~ Sex + age + diagnosis_summarized + time_from_ED_dich + lines_pretreatment_cat+ material +tumor_infiltration +Treatment_type,  data=df_invivo, overall="Total",  render.continuous=my.render.cont, render.categorical=my.render.cat)
Total
(N=80)
Sex
Male 48 (60 %)
Female 32 (40 %)
Age (years)
Median [Min, Max] 68.5 [20.0, 91.0]
Diagnosis
AML 34 (42 %)
ALL 2 (2 %)
DLBCL 4 (5 %)
Burkitt lymphoma 1 (1 %)
Follicular lymphoma 3 (4 %)
Mantle cell lymphoma 5 (6 %)
CLL 25 (31 %)
B-PLL 1 (1 %)
T-PLL 4 (5 %)
T-cell lymphoma, NOS 1 (1 %)
Time from initial diagnosis (months)
≤12 49 (61 %)
>12 31 (39 %)
Prior lines of treatment
≥3 6 (8 %)
0 54 (68 %)
1 11 (14 %)
2 9 (11 %)
Tumor sample origin
Periphal blood 59 (74 %)
Bone marrow 12 (15 %)
Lymph node 9 (11 %)
Tumor infiltration of sample
Median [Min, Max] 84.5 [52.0, 99.0]
Prescribed treatment at study entry
Chemotherapy 28 (35 %)
Chemotherapy + small molecule inhibitors 7 (9 %)
Immuno-Chemotherapy 16 (20 %)
Immunotherapy 2 (2 %)
Small molecule inhibitors 27 (34 %)

4.2 4.2 Patient list

For appendix: includes clinical data for each patient

df_patients_supp<-df_invivo%>%select(patientID, diagnosis, age, Sex, material, tumor_infiltration, lines_pretreatment, treatment_spec, Response_def_parameter_final,  Resp_protocol, Death, time_EFS)
df_patients_supp<-df_patients_supp%>%mutate(time_EFS=ifelse(Resp_protocol=="na", "na", time_EFS))

df_patients_supp %>% DT::datatable(width = 700)
# write_xlsx(df_patients_supp,"04_Figures/table_S3_patients.xlsx")

remove(df_patients_supp)

4.3 4.3 Median Follow-up

Calculated according to the reverse Kaplan-Meier method

fit_FU <- survfit(Surv(time_OS/365, Death_r) ~ 1, data=df_invivo)
surv_median(fit_FU)
## Warning: `select_()` was deprecated in dplyr 0.7.0.
## Please use `select()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
##   strata median lower upper
## 1    All   2.11  1.95  2.62
remove(fit_FU)

5 5 Primary endpoint

5.1 5.1 Calculation of primary endpoint

df_invivo%>%summarize(
         mean.timelabreport=mean(time_report),
         median.timelabreport=median(time_report),
         min.timelabreport=min(time_report),
         max.timelabreport=max(time_report))
## # A tibble: 1 x 4
##   mean.timelabreport median.timelabreport min.timelabreport max.timelabreport
##   <drtn>             <drtn>               <drtn>            <drtn>           
## 1 3.98 days          3 days               2 days            17 days
#calculate rate with labreport during 7 days
df_prim_endpoint<-df_invivo%>% select(patientID, time_report) %>% dplyr::mutate(within7d=ifelse(time_report<=7, 1, 0)) %>% dplyr::mutate(time_to_analysis=ifelse(time_report<=7, "within 7 days", "more than 7 days")) %>% 
   dplyr::mutate(duration_analysis_cat=case_when(time_report<=3 ~ "within 3 days", (time_report>4 &time_report<=7) ~ "within 4-7 days", (time_report>=8 &time_report<=10) ~ "within 8-10 days", time_report>10 ~ "more than 10 days" ))

bintest <- binom.test(sum(df_prim_endpoint$within7d), nrow(df_prim_endpoint))
bintest
## 
##  Exact binomial test
## 
## data:  sum(df_prim_endpoint$within7d) and nrow(df_prim_endpoint)
## number of successes = 73, number of trials = 80, p-value = 6e-15
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.828 0.964
## sample estimates:
## probability of success 
##                  0.912
cat("For ", sum(df_prim_endpoint$within7d)," of ", length(df_prim_endpoint$within7d), "patients the primary endpoint was reached. Therefore, for ", (bintest$estimate*100),"% (",round(bintest$conf.int[1],3),",", round(bintest$conf.int[2],3),") of participants ex-vivo drug testing was performed within a week.")
## For  73  of  80 patients the primary endpoint was reached. Therefore, for  91.2 % ( 0.828 , 0.964 ) of participants ex-vivo drug testing was performed within a week.
remove(bintest)

#calculate rate with labreport during 3 days
df_invivo%>% select(patientID,time_report) %>% dplyr::mutate(within3d=ifelse(time_report<4, 1, 0)) %>% dplyr::count(within3d) %>% mutate(percent = n / sum(n)*100)
## # A tibble: 2 x 3
##   within3d     n percent
##      <dbl> <int>   <dbl>
## 1        0    30    37.5
## 2        1    50    62.5

5.2 5.2 Visualisation primary endpoint (barplot)

df_prim_endpoint<-df_prim_endpoint%>% dplyr::count(time_to_analysis) %>% mutate(percent = n / sum(n)*100)

df_prim_endpoint$time_to_analysis<-factor(df_prim_endpoint$time_to_analysis, levels=c("within 7 days", "more than 7 days"), labels=c("within 7 days", "more than\n7 days"))

ggplot(data=df_prim_endpoint, aes(x=as.factor(time_to_analysis), y=percent, fill=time_to_analysis))+
  geom_bar(stat="identity", width=0.5, alpha=0.85)+
  xlab("")+
  annotate("text", x=1, y=98, label= "73/80", size=6)+
  annotate("text", x=2, y=15, label= "7/80", size=6)+
  scale_y_continuous("Number of patients in %", limits = c(0, 100), breaks=c(0, 20, 40, 60, 80, 100), expand = expansion(mult = c(0, 0.05)))+
  scale_fill_manual(values=c(red1, blue1), name="Time to report")+
  theme_figures_wo_legend

# ggsave(width = 6, height = 4, path=(paste0(filepath,"04_Figures")), filename = "Figure_2B_Time_to_report.svg")


# remove(df_prim_endpoint)

6 6 Pathways and drug difference between in-vivo response groups in chemotherapy patients

6.1 6.1 Prepare base data frame

df_chemo$Pathway <-factor(df_chemo$Pathway, 
                                levels=c("ABL (BCR)", "ALK", "Apoptosis", "BCR", "Bromodomain", "CDK", "Cell cycle" , 
                                         "Chemotherapy", "DDR", "FLT3", "Hedgehog", "HGF", "Histone acetyltransferase", 
                                         "Histone deacetylase", "Histone methyltransferase", "IDH", "JAK/STAT", "MAPK", 
                                         "Notch", "Nuclear export", "PI3K/AKT/mTOR", "PKC", "PLK", "Promiscuous" , "Proteasome", 
                                         "Stress response", "TLR", "TNF/NFKB" ), 
                                labels=c("BCR/ABL", "ALK", "Apoptosis", "B cell receptor", "Bromodomain", "Cyclin dependent kinases", 
                                         "Cell cycle control" , "Chemotherapy", "DNA damage response", "FLT3", "Hedgehog", "HGF", 
                                         "Histone acetyltransferase", "Histone deacetylase", "Histone methyltransferase", "IDH", 
                                         "JAK/STAT", "MAPK", "Notch", "Nuclear export", "PI3K/AKT/mTOR", "PKC", "PLK", "Other" , 
                                         "Proteasome", "Stress response", "TLR", "TNF/NFKB" ))

df_chemo_bckup <- rename(df_chemo, response = resp_protocol)

df_chemo <- filter(df_chemo, resp_protocol != "SD") %>%
  #rename some cols
  rename(response = resp_protocol) %>%
  #AUC per pathway per patient
  group_by(patientID, Pathway) %>% mutate(AUC_pathway_patient = mean(AUC)) %>%
  ungroup() 

6.2 6.2 Perform statistical tests and figures

6.2.1 6.2.1 t-tests by drug, R vs PD

  #eliminate conc-specific rows/cols          
df_p_drugAUC <- select(df_chemo, 
                       -concentration, -concIndex, -normVal, -AUC_pathway_patient) %>% unique() %>%
  #"AUC" to "metric" for functions
  mutate(metric = AUC) %>% 
  #perform standard two-sample t-test
  group_by(name) %>% nest() %>%
  mutate(pval = map(data, t_test) ) %>% 
  unnest( c(data, pval) ) %>% ungroup() %>%
  #how many drugs per pathway
  group_by(Pathway) %>% mutate( n_drugs_pathway = length(unique(name)) ) %>% 
  #How many significant drugs per pathway (raw p 0.05)
  mutate( n_sign_drugs_pathway = length( unique(name[pval < 0.05]) ) ) %>%
  #Prop of significant drugs per pathway
  mutate(prop_sig_pathway = (n_sign_drugs_pathway/n_drugs_pathway)) %>%
  ungroup() %>%
  #arrange by pathway with highest prop and by most significant drugs
  arrange( desc(prop_sig_pathway), Pathway, pval )

df_p_drugAUC %>% select(Pathway, name, prop_sig_pathway, pval) %>% unique() %>% 
  mutate_if(is.numeric, formatC, digits=2) %>% DT::datatable(width = 700)

6.3 6.3 Figure 3A: Association between pathway ex vivo response and in vivo response in patients treated with chemotherapy.

6.3.1 6.3.1 Prepare data

#mean per pathway per in-vivo group
plotTab <- group_by(df_p_drugAUC, Pathway, response) %>% 
  summarise(AUC_pathway_response = mean(AUC), .groups="keep") %>% 
  #difference R vs PD
  group_by(Pathway) %>%
  summarise(diff_RvsPD = AUC_pathway_response[response == "R"] - AUC_pathway_response[response == "PD"], 
            direction = ifelse(diff_RvsPD < 0,"decreased","increased"), .groups = "keep") %>% ungroup() %>%
  #arrange and update arrangement for plot
  arrange( diff_RvsPD ) %>% mutate(Pathway = factor(Pathway, levels = unique(Pathway)))

Pathway_order<-plotTab %>% 
  select(Pathway) %>% 
  mutate(Pathway=as.character(Pathway))%>% 
  unlist()

6.3.2 6.3.2 Plot waterfall

plot_3A<-ggplot(plotTab, aes(x=Pathway, y=diff_RvsPD, fill=direction)) +
  geom_col(color="black", size=0.7, alpha=0.65, ) +
  scale_fill_manual(values = c("decreased"=red1,"increased"=blue1), name = "Viability")+ 
  coord_cartesian(ylim = c( min((plotTab$diff_RvsPD)) , 0.2 ) ) +
  ylab("Viability difference\n") + xlab("") +
  theme_figures_x_angle_45_wo_legend+
  ggtitle("\u0394 = viability (R) – viability (PD)")


plot_3A

# remove(plotTab)

6.4 6.4 Panel B: Significance in R vs PD comparisons by drug (t-tests).

6.4.1 6.4.1 Prepare data

plotTab <- select(df_p_drugAUC, Pathway, name, pval) %>% unique() %>% 
  mutate(log10_pval = -log10(pval), significant = ifelse(pval < 0.05,T,F),
         significant = factor(significant, levels=unique(significant) ),
         name = as.character(name)) 

plotTab[plotTab$name == "Abemaciclib (LY2835219)",]$name <- "Abemaciclib"

# drug_order<-plotTab%>% 
#   arrange(pval)%>%
#   select(Pathway)%>%
#   mutate(Pathway=as.character(Pathway))%>% 
#   unlist()%>% 
#   unique()

6.4.2 6.4.2 Plot

plot_3B<-ggplot(plotTab, aes(x=factor(Pathway, levels = Pathway_order), y=log10_pval, fill=significant)) + 
  geom_point(shape=21, color ="black", stroke=0.5, size=2.5, alpha=0.6, show.legend = F) +
  scale_fill_manual(values = c("TRUE"=red1, "FALSE"="white")) +
  geom_hline(yintercept = -log10(0.05), linetype="22", col="gray37") +
  gghighlight(significant == T, use_direct_label = F, 
              unhighlighted_params = list(fill = NULL, color="gray") ) +
  geom_text_repel(data = filter(plotTab, significant == T), aes(label=name), 
                   size=fontsize_small-4, max.overlaps = Inf, min.segment.length = 0.4, 
                  box.padding = 0.65, segment.colour="grey")+
  ylab(expression(-log[10]*'('*italic(P)~value*')')) +
  xlab("") +
 theme_classic()+
  theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
        axis.text.x =element_text(size=fontsize, angle=45, vjust=1, hjust=1),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        panel.grid.major.x = element_line(size = 0.5,linetype = "13"), 
        axis.ticks.x = element_line(size = 0.5, color = "black", linetype = "13"), 
        axis.ticks.length.x = unit(0.3, "cm"),
        legend.position = "none")
 
plot_3B 

# remove(plotTab)

6.5 6.5 Panel C: Dose-response curves and AUC box-plots for exemplary drugs.

Selected drugs

example_drugs <- c("Idarubicin", "Vindesine", "Palbociclib", "Ganetespib")

Fit curves

#data storage
df_smooth_curve <- data.frame()

for ( drugname in example_drugs ){
  #prepare data
  df <- filter(df_chemo, name == drugname ) %>% 
    select(patientID, response, name, concentration, normVal) %>%
    arrange( desc(response) ) %>% mutate(response = factor(response, levels = unique(response)))
  #call fun
  res <- suppressWarnings(fit_log_curve(df, drugname))
  #store data (1k predicted viability values for plotting a smooth curve)
  df_smooth_curve <- rbind(df_smooth_curve, res[[5]] )
}

prepare data

plotTab <- select(df_chemo, Pathway, name, patientID, response, concentration, normVal) %>%
  #PD points on top
  arrange( desc(response) ) %>% mutate(response = factor(response, levels=unique(response)))

make figs

plotList_DResp <- lapply( example_drugs ,function(drugname){
  #plot values: real with dots, predicted with line
  ggplot(NULL,aes(color=response, fill=response)) + 
    geom_line(data = df_smooth_curve[df_smooth_curve$name == drugname,], 
              aes(x= concentration, y= viab_pred), size=1) + 
    geom_point(data = plotTab[plotTab$name == drugname,] , 
               aes(x= concentration, y=normVal), 
               shape = 21, size=3.5, alpha = 0.5, color="black") +
    scale_color_manual(values = c("R" = blue1, "PD"=red2)) +
    scale_fill_manual(values = c("R" = blue2, "PD"=red1)) +
    
    coord_cartesian(ylim = c(0,1.5)) +
    scale_x_continuous(trans='log10') +
    
    labs(title =  paste0("",drugname) )+
    xlab(expression("Concentration"~(mu*M))) + ylab("Viability") +
    guides(fill="none", color="none") +
    theme_classic()+
    theme(plot.title=element_text(size=fontsize+4, hjust = 0, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
        axis.text.x =element_text(size=fontsize),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        axis.line = element_blank(), 
        panel.border = element_rect(size=1.5, fill = NA),
        axis.ticks = element_line(size=1.5),
        legend.text = element_text(size=fontsize),
        legend.title = element_text(size=fontsize+2))
})


plotList_box <- lapply( example_drugs ,function(drugname){
  #get the df
  df_plot <- filter(df_p_drugAUC, name == drugname)
  #make plot
  p <- ggplot(df_plot, aes(x=response, y=AUC)) + 
    geom_boxplot(aes(fill = response), width = 0.8, size=0.8, alpha = 0.85, 
                 color = "black", outlier.shape = NA) +
    scale_fill_manual(values = c("PD"=red1,"R"=blue1) ) +
    geom_beeswarm(aes(color = response),cex=3, dodge.width=0.8, shape=21, 
                  fill="slategray4", size=4, alpha= 0.7) +
    scale_color_manual(values = c("black", "black") ) + 
    geom_hline(yintercept = 1, linetype="22", col="gray30", size=0.7) +
    coord_cartesian(ylim = c(0,1.5)) +
    
    labs(title= "" ) +
    xlab(expression(paste(italic("in vivo")," response"))) +
    ylab("Viability (AUC)") +
    guides(color="none", fill="none") +
    theme_classic()+
    theme(plot.title=element_text(size=fontsize+4, hjust = 0, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
        axis.text.x =element_text(size=fontsize),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        axis.line = element_blank(), 
        panel.border = element_rect(size=1.5, fill = NA),
        axis.ticks = element_line(size=1.5),
        legend.text = element_text(size=fontsize),
        legend.title = element_text(size=fontsize+2))

  return( p )
})

#merge both plot lists
plotList <- c(rbind(plotList_DResp, plotList_box))

get legend

plot figs

design<-"
A#B
C#D
E#F
G#H
III
"

plot_3C<-
plotList[[1]]+
plotList[[2]]+
plotList[[3]]+
plotList[[4]]+
plotList[[5]]+
plotList[[6]]+
plotList[[7]]+
plotList[[8]]+
  legend_3C+
  
  plot_layout(design = design, widths = c(1,0.02,0.98))

plot_3C

remove(df, df_smooth_curve, plotList, plotList_box, plotList_DResp, plotTab, res, legend_3C)

7 7 Elastic net regression model

# Set a seed for reproducibility
set.seed(20191120)


#Modify drug names for formula function
tmp <- df_drugs %>% mutate(Drug = as.character(Drug), Drug = case_when(Drug == "8-Azaguanin" ~ "Azaguanin", Drug == "Azacytidine (Vidaza)" ~ "Vidaza",Drug == "Decitabine (Dacogen)" ~ "Decitabine", TRUE ~ Drug))
chemotherapies <- (tmp %>% filter(pathway == "Chemotherapy", !(Drug %in% c("Cisplatin", "Cyclophasphamide"))))$Drug

lasso_analysis_set <- data_all %>%
  # filter only chemopatients
  filter(chemo_pat==TRUE)%>%
  filter(!is.na(Treatment_type), Treatment_type!="na", Resp_protocol %in% c("R", "PD"))%>%
  # filter only patients with an evaluable in vivo response
  filter(Response_not_evaluable!="1")%>%
  # filter patients only with both plates available
  filter(Count_plates==2) %>%
  #restrict to chemotherapies, modify drug names
  mutate(Drug = as.character(Drug)) %>%
  mutate(Drug = case_when(Drug == "8-Azaguanin" ~ "Azaguanin", Drug == "Azacytidine (Vidaza)" ~ "Vidaza",Drug == "Decitabine (Dacogen)" ~ "Decitabine", TRUE ~ Drug)) %>%
  filter(Drug %in% chemotherapies)  %>%
  select(patientID, Drug, AUC, Resp_protocol)

#Transform data to wide format for regression analysis
lasso_analysis_set <- lasso_analysis_set %>% tidyr::spread("Drug", AUC)

formula <- as.formula(paste("Resp_protocol ~ ", paste(chemotherapies, collapse = " + "), sep = ""))

x.bestr <- model.matrix(formula, lasso_analysis_set)
y.bestr <- lasso_analysis_set$Resp_protocol

lasso.bestr.models <- repeated_glmnet(X = x.bestr, Y = y.bestr, alpha = 0.3, nlambda = 100, nfolds = 3)

list.results.bestr <- give_repeated_lasso_results(lasso.bestr.models, type.measure ="auc")

# 
# print.xtable(xtable(list.results.bestr[[1]]), include.rownames=FALSE, table.placement = "ht", type = "html")
# write.csv(list.results.bestr[[1]], file = "elastic_net_results.csv")
knitr::kable(list.results.bestr[[1]][,-1])
medianOR Proportion_selected median_AUROC
Vincristine 0.959 1.000 0.847
Idarubicin 0.967 0.987 0.847
Vindesine 0.971 0.985 0.847
Mitoxantrone 0.974 0.969 0.847
Cladribine 0.973 0.884 0.844
Fludarabine 0.976 0.532 0.833
Azaguanin 0.965 0.445 0.836
Etoposide 0.973 0.445 0.836
Thioguanin 0.933 0.389 0.839
Pentostatin 1.120 0.262 0.844
Bendamustine 0.941 0.249 0.847
Pralatrexate 0.921 0.249 0.847
Nelarabine 1.009 0.150 0.861
Clofarabine 1.141 0.145 0.865
Daunorubicin 1.128 0.138 0.867
Cytarabine 1.074 0.122 0.867
Doxorubicin 1.031 0.112 0.867
Chlorambucil 1.010 0.107 0.876
Vidaza 1.012 0.101 0.867
Decitabine 1.069 0.014 0.877
(Intercept) 1.000 0.000 0.500
tabforFig3D<-list.results.bestr[[1]][1:6,-1]
colnames(tabforFig3D)<-c("Median OR", "Selected Proportion", "Median AUROC")

corrplot(list.results.bestr[[2]], type="upper", col=brewer.pal(n=8, name="RdYlBu"))

tabforFig3D<-tabforFig3D %>% 
  mutate(`Median OR`=round(`Median OR`,2),
         `Selected Proportion`=round(`Selected Proportion`,2),
         `Median AUROC`=round(`Median AUROC`,2))

Fig3D<-  ggplot(rownames_to_column(tabforFig3D, var = "Drug"))+
   annotate(geom = "table", size=fontsize_small-4, x = 0,y = 0, label = list(cbind(rownames_to_column(tabforFig3D, var = "Drug"))))+
    theme_void()

Fig3D

8 8 EFS analyses

8.1 8.1 Drugs which were relevant in elastic net model

tmp<-data_all%>%
  # filter only chemopatients
  filter(chemo_pat==TRUE)%>%
  filter(!is.na(Treatment_type), Treatment_type!="na")%>%
  # filter only patients with an evaluable in vivo response
  filter(Response_not_evaluable!="1")%>%
  # filter only relevant drugs from elastic net
  filter(Pathway=="Chemotherapy")%>%filter(Drug=="Vindesine" | Drug=="Idarubicin" |  Drug=="Mitoxantrone" |Drug=="Vincristine" |Drug=="Cladribine"|Drug=="Fludarabine")%>%
  # filter patients only with both plates available (to be consistent with elastic net model)
  filter(Count_plates==2)

tmp$Drug<- as.factor(tmp$Drug)

tmp<-tmp%>%select(patientID, Drug, AUC, Event, time_EFS)

# AUC is scaled such that a unit change of the regressor corresponds to 10% change in cell viability. 
tmp<-tmp%>%mutate(AUC_percent=AUC*100)
tmp<-tmp%>%mutate(AUC_10percent=AUC*100/10)

# loop to calculate HR + 95%CI + p for each drug
result<-data.frame(matrix(nrow =1, ncol=5))
colnames(result)<- c("drug", "HR",  "lower", "upper","p_value_cox")

for(i in unique(tmp$Drug)){
tmp2<-tmp%>%filter(Drug==i)
cox<-coxph(Surv(time_EFS/365, Event) ~ (AUC_10percent), 
               data=tmp2)
cox_sum<-summary(cox)


result[i, 1] <-i
result[i, 2] <-cox_sum[["coefficients"]][, "exp(coef)"]
result[i, 3] <-cox_sum[["conf.int"]][, "lower .95"]
result[i, 4] <-cox_sum[["conf.int"]][, "upper .95"]
result[i, 5] <-cox_sum[["coefficients"]][, "Pr(>|z|)"]
}

result<-result[-which(rownames(result)==1),]
tablefig3D<-result[-which(rownames(result)==1),]

# plot results
# knitr::kable(tablefig3D)

remove(tmp, tmp2, cox, cox_sum)

Plotting of Forest plot

result$drug<-factor(result$drug, levels=rev(c("Vincristine", "Idarubicin","Vindesine", "Mitoxantrone", "Cladribine", "Fludarabine")))

plot_3E<-ggplot(data=result ,aes(x=drug, y=HR, ymin=lower, ymax=upper,
                     )) +  geom_pointrange() + 
     # scale_x_discrete(limits=result$drug,breaks=rev(levels(result$drug)), name="")+
     scale_y_log10() +
  xlab("")+
     scale_color_manual(values="black", name="")+
     coord_flip() + 
     geom_hline(aes(yintercept=1), colour="black", size=1.5,
                linetype="dashed", alpha=0.3)  +
     # annotate("text", x = 1:nrow(result)+0.2, y = result$HR+0.003,
     #         label = ifelse( result$p_value_cox<0.001, paste0("p<","0.001"), 
     #            paste0("p=", round(result$p_value_cox,3) ) ), size=fontsize_small-4)+
  geom_text(nudge_x =.3, label=paste0("p=", round(result$p_value_cox,3)), size=fontsize_small-5)+
    theme(panel.grid.minor = element_blank(),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
        axis.text.x =element_text(size=fontsize),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        legend.position = "none",
        panel.background = element_rect(fill = "white", color="black"),
        panel.grid.major = element_blank(),
        axis.ticks.x = element_line(size=1.5),
        axis.ticks.y = element_blank())+
  ylab("Hazard Ratio")

plot_3E  

# remove(result)

8.2 8.2 EFS-KM curves for selected drugs

8.2.1 8.2.1 Vindesine ~ EFS in chemogroup

x<-"Vindesine"

# prepare data
tmp<-data_all%>%
  # filter only chemopatients
  filter(chemo_pat==TRUE)%>%
  filter(!is.na(Treatment_type), Treatment_type!="na")%>%
  # filter only patients with an evaluable in vivo response
  filter(Response_not_evaluable!="1")%>%
  # filter only relevant drugs from elastic net
  filter(Pathway=="Chemotherapy")%>%filter(Drug=="Vindesine" | Drug=="Idarubicin" |  Drug=="Mitoxantrone" |Drug=="Vincristine" |Drug=="Cladribine"|Drug=="Fludarabine")%>%
  # filter patients only with both plates available (to be consistent with elastic net model)
  filter(Count_plates==2)%>%
  filter(Drug==x)

# perform maximally selected log rank statistics
m_EFS<-maxstat.test(Surv(time_EFS/365, Event) ~ AUC, data=tmp,
    smethod="LogRank", pmethod="Lau94", minprop = 0.2, 
 
                                maxprop = 0.8)
m_EFS
## 
## Maximally selected LogRank statistics using Lau94
## 
## data:  Surv(time_EFS/365, Event) by AUC
## M = 2, p-value = 0.1
## sample estimates:
## estimated cutpoint 
##              0.766
# divide patients into two groups based on cutpoint
tmp$drug_sensitive<-tmp$AUC<=m_EFS$estimate

#  perform survival analyses
fit_EFS <- survfit(Surv(time_EFS/365, Event) ~ drug_sensitive, data=tmp)
surv_median(fit_EFS)
##                 strata median  lower upper
## 1 drug_sensitive=FALSE 0.0932 0.0521 0.392
## 2  drug_sensitive=TRUE 0.4849 0.2301    NA
survdiff(Surv(time_EFS/365, Event) ~ drug_sensitive, data=tmp)
## Call:
## survdiff(formula = Surv(time_EFS/365, Event) ~ drug_sensitive, 
##     data = tmp)
## 
##                       N Observed Expected (O-E)^2/E (O-E)^2/V
## drug_sensitive=FALSE 15       14     7.29      6.19      8.44
## drug_sensitive=TRUE  28       19    25.71      1.75      8.44
## 
##  Chisq= 8.4  on 1 degrees of freedom, p= 0.004
# plot KM curve
ggsurvplot_EFS_Vind<-ggsurvplot(fit_EFS,  conf.int =FALSE, risk.table = TRUE, pval = FALSE,
           title="Vindesine",
           xlab ="Time (in years)",
           xlim = c(0, 2.5),
          break.time.by=0.5,
           ylab ="\nEvent-free survival probability\n",
           surv.median.line = "hv",
           tables.theme=clean_theme(),
           legend="none",
           legend.title=(""),
           legend.labs=c("less sensitive", "sensitive"),
           palette = c(blue1, red1),
           ggtheme=theme_figures)

ggsurvplot_EFS_Vind

plot_EFS_Vind<-ggsurvplot_EFS_Vind$plot

remove(tmp, x, m_EFS, fit_EFS, ggsurvplot_EFS_Vind)

8.2.2 8.2.2 Vincristine ~ EFS in chemogroup

x<-"Vincristine"

# prepare data
tmp<-data_all%>%
  # filter only chemopatients
  filter(chemo_pat==TRUE)%>%
  filter(!is.na(Treatment_type), Treatment_type!="na")%>%
  # filter only patients with an evaluable in vivo response
  filter(Response_not_evaluable!="1")%>%
  # filter only relevant drugs from elastic net
  filter(Pathway=="Chemotherapy")%>%filter(Drug=="Vindesine" | Drug=="Idarubicin" |  Drug=="Mitoxantrone" |Drug=="Vincristine" |Drug=="Cladribine"|Drug=="Fludarabine")%>%
  # filter patients only with both plates available (to be consistent with elastic net model)
  filter(Count_plates==2)%>%
  filter(Drug==x)

# perform maximally selected log rank statistics
m_EFS<-maxstat.test(Surv(time_EFS/365, Event) ~ AUC, data=tmp,
    smethod="LogRank", pmethod="Lau94", minprop = 0.2, 
                             maxprop = 0.8)
m_EFS
## 
## Maximally selected LogRank statistics using Lau94
## 
## data:  Surv(time_EFS/365, Event) by AUC
## M = 2, p-value = 0.3
## sample estimates:
## estimated cutpoint 
##              0.722
# divide patients into two groups based on cutpoint
tmp$drug_sensitive<-tmp$AUC<=m_EFS$estimate

#  perform survival analyses
fit_EFS <- survfit(Surv(time_EFS/365, Event) ~ drug_sensitive, data=tmp)
surv_median(fit_EFS)
##                 strata median  lower upper
## 1 drug_sensitive=FALSE 0.0932 0.0767 0.603
## 2  drug_sensitive=TRUE 0.4849 0.2082    NA
survdiff(Surv(time_EFS/365, Event) ~ drug_sensitive, data=tmp)
## Call:
## survdiff(formula = Surv(time_EFS/365, Event) ~ drug_sensitive, 
##     data = tmp)
## 
##                       N Observed Expected (O-E)^2/E (O-E)^2/V
## drug_sensitive=FALSE 14       13     7.22      4.64      6.29
## drug_sensitive=TRUE  29       20    25.78      1.30      6.29
## 
##  Chisq= 6.3  on 1 degrees of freedom, p= 0.01
# plot KM curve
ggsurvplot_EFS_Vinc<-ggsurvplot(fit_EFS,  conf.int =FALSE, risk.table = TRUE, pval = FALSE,
           title="Vincristine",
           xlab ="Time (in years)",
           xlim = c(0, 2.5),
           break.time.by=0.5,
           ylab ="\nEvent-free survival probability\n",
           surv.median.line = "hv",
           tables.theme=clean_theme(),
           legend="none",
           legend.title=(""),
           legend.labs=c("less sensitive", "sensitive"),
           palette = c( blue1, red1),
           ggtheme=theme_figures)

ggsurvplot_EFS_Vinc

plot_EFS_Vinc<-ggsurvplot_EFS_Vinc$plot

remove(tmp, x, m_EFS, fit_EFS, ggsurvplot_EFS_Vinc)

8.2.3 8.2.3 Cladribine ~ EFS in chemogroup

x<-"Cladribine"

# prepare data
tmp<-data_all%>%
  # filter only chemopatients
  filter(chemo_pat==TRUE)%>%
  filter(!is.na(Treatment_type), Treatment_type!="na")%>%
  # filter only patients with an evaluable in vivo response
  filter(Response_not_evaluable!="1")%>%
  # filter only relevant drugs from elastic net
  filter(Pathway=="Chemotherapy")%>%filter(Drug=="Vindesine" | Drug=="Idarubicin" |  Drug=="Mitoxantrone" |Drug=="Vincristine" |Drug=="Cladribine"|Drug=="Fludarabine")%>%
  # filter patients only with both plates available (to be consistent with elastic net model)
  filter(Count_plates==2)%>%
  filter(Drug==x)

# perform maximally selected log rank statistics
m_EFS<-maxstat.test(Surv(time_EFS/365, Event) ~ AUC, data=tmp,
    smethod="LogRank", pmethod="Lau94", minprop = 0.2, 
                             maxprop = 0.8)
m_EFS
## 
## Maximally selected LogRank statistics using Lau94
## 
## data:  Surv(time_EFS/365, Event) by AUC
## M = 2, p-value = 0.5
## sample estimates:
## estimated cutpoint 
##              0.548
# divide patients into two groups based on cutpoint
tmp$drug_sensitive<-tmp$AUC<=m_EFS$estimate

#  perform survival analyses
fit_EFS <- survfit(Surv(time_EFS/365, Event) ~ drug_sensitive, data=tmp)
surv_median(fit_EFS)
##                 strata median  lower upper
## 1 drug_sensitive=FALSE  0.132 0.0438    NA
## 2  drug_sensitive=TRUE  0.353 0.0986  1.66
survdiff(Surv(time_EFS/365, Event) ~ drug_sensitive, data=tmp)
## Call:
## survdiff(formula = Surv(time_EFS/365, Event) ~ drug_sensitive, 
##     data = tmp)
## 
##                       N Observed Expected (O-E)^2/E (O-E)^2/V
## drug_sensitive=FALSE  9        9     4.68     3.982      4.92
## drug_sensitive=TRUE  34       24    28.32     0.658      4.92
## 
##  Chisq= 4.9  on 1 degrees of freedom, p= 0.03
# plot KM curve
ggsurvplot_EFS_Clad<-ggsurvplot(fit_EFS,  conf.int =FALSE, risk.table = TRUE, pval = FALSE,
           title="Cladribine",
           xlab ="Time (in years)",
           xlim = c(0, 2.5),
           break.time.by=0.5,
           ylab ="\nEvent-free survival probability\n",
           surv.median.line = "hv",
           tables.theme=clean_theme(),
           legend="none",
           legend.title=(""),
           legend.labs=c("less sensitive", "sensitive"),
           palette = c( blue1, red1),
           ggtheme=theme_figures)

ggsurvplot_EFS_Clad

plot_EFS_Clad<-ggsurvplot_EFS_Clad$plot

remove(tmp, x, m_EFS, fit_EFS, ggsurvplot_EFS_Clad)

8.2.4 8.2.4 Build plot with KM curves of selected drugs (Figure 3F)

design1<-"
AB
DD
"

plot_3F<-
  plot_EFS_Vinc+
  plot_EFS_Vind+
  # plot_EFS_Clad+
  legend_3F+
  plot_layout(design=design1, heights = c(1,0.15))
plot_3F

# remove(plot_EFS_Clad, plot_EFS_Vinc, plot_EFS_Vind, legend_3F)

9 9 Plot Figure 3

library("Cairo")


design<-"
AAAC
DDEE
FFFF
"
tp<-theme(plot.tag = element_text(size=22))

  wrap_elements(plot_3A+ plot_3B+ plot_layout(ncol=1)+ plot_annotation(tag_levels = list(c('A', 'B'))))+
  wrap_elements(plot_3C+ plot_annotation(tag_levels = list(c('C'))))+
  wrap_elements( Fig3D +tp + plot_annotation(tag_levels = list(c('D'))))+
  wrap_elements(plot_3E+ plot_annotation(tag_levels = list(c('E'))))+
  wrap_elements(plot_3F+ plot_annotation(tag_levels = list(c('F'))))+ 
  plot_layout(design=design, widths = c(1,1,.2,1.8), heights = c(2,.5,.7))+
    plot_annotation(title="Figure 3") & theme(title = element_text(size=30))

# ggsave(path=(paste0(filepath,"04_Figures")), filename = "Figure_3.svg")
  
# ggsave(path=(paste0(filepath,"04_Figures")), filename = "Figure_3.pdf", device = cairo_pdf)

10 10 Burkitt lymphoma patient example (S005)

Burkitt Lymphoma: pretreatment with GMALL, during study treatment with R-DHAP –> PD with new liver lesions, then salvage treatment with MTX –> CR

# Prepare data: filter only for patient S005, label specific drugs which the patient had received to use these for colorcoding later in the plot
df_S005<-data_all%>%filter(patientID=="S005")%>%mutate(Drug_coding=case_when((Drug=="Doxorubicin" | Drug=="Cytarabine" ) ~ "refractory to treatment", Drug=="Pralatrexate" ~ "response to salvage treatment", (Drug!="Doxorubicin" &Drug!="Cytarabine" & Drug!="Pralatrexate") ~ "treatment not administered"))

# rename Drugs which includes the brand name in the original drug list
df_S005<-df_S005%>%mutate(Drug=ifelse(Drug=="Azacytidine (Vidaza)", "Azacytidine", Drug))
df_S005<-df_S005%>%mutate(Drug=ifelse(Drug=="Decitabine (Dacogen)", "Decitabine", Drug))

# define theme for figure
theme_figures_S005<- theme_classic()+
  theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
         axis.text.x =element_text(size=fontsize, angle=45, vjust=1, hjust=1),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        axis.line = element_line(size=0.6),
        axis.ticks = element_line(size=1.5),
        legend.text = element_text(size=fontsize),
        legend.title = element_text(size=fontsize+2),
        legend.position ="bottom")


plot_S005_ESS<-df_S005%>%filter(Pathway=="Chemotherapy")%>%
  ggplot(aes(x=reorder(Drug, ESS_AUC), y=ESS_AUC, fill=Drug_coding)) +
  geom_bar(stat="identity", alpha=0.85)+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Effect size score\n", limits=c(-0.5, 0.75))+
  geom_hline(yintercept = 0)+
  scale_fill_manual(name=expression(paste(italic("in vivo")," response")), values=c(blue1, red1, "gray84"))+
  theme_figures_S005+
  theme(axis.text.x =element_text(angle=90, vjust=0.25))+
  guides(fill = guide_legend(nrow = 3))      

plot_S005_ESS

# ggsave( path=(paste0(filepath,"04_Figures")), filename = "Figure_4A_Burkitt_lymphoma.svg")

remove(df_S005, plot_S005_ESS, theme_figures_S005)

11 11 CLL patients with venetoclax and ibrutinib treatment

11.1 11.1 Ibrutinib

theme_CLL_ibr<- theme_classic()+
  theme(plot.title=element_text(size=fontsize+4, hjust = 0.5, face="bold"),
        text=element_text(size=fontsize),
        axis.text.y =element_text(size=fontsize),
         axis.text.x =element_text(size=fontsize, angle=45, vjust=1, hjust=1),
        axis.title.y =element_text(size=fontsize+2),
        axis.title.x =element_text(size=fontsize+2),
        axis.line = element_line(size=0.6),
        axis.ticks = element_line(size=1.5),
        legend.text = element_text(size=fontsize),
        legend.title = element_text(size=fontsize+2),
        legend.position = c(0.8, 0.2))

# prepare CLL ibrutinib dataset
df_CLL_ibr<-data_all%>%mutate(Ibrutinib=ifelse(grepl("Ibrutinib", treatment_spec), 1, 0))
df_CLL_ibr<-df_CLL_ibr%>%filter(Drug=="Ibrutinib")%>%filter(Ibrutinib==1)%>%filter(diagnosis=="CLL")%>%filter(Response_not_evaluable==0)
df_CLL_ibr$Resp_protocol <-factor(df_CLL_ibr$Resp_protocol, levels=c("R", "SD"), 
                       labels=c("response","stable disease")) # patients had only R or SD, no patients with PD


plot_ibr<-df_CLL_ibr%>%ggplot(aes(x=reorder(patientID, -ES_AUC), y=-ES_AUC, fill=Resp_protocol)) +
  geom_bar(stat="identity", alpha=0.75)+
  scale_x_discrete(name="\nCLL patients")+
  scale_y_continuous(name=expression(paste("Effect size (",italic("ex vivo"),")")), limits=c(-0.4, 0), breaks=c( -0.4, -0.2, 0), labels=c( 0.4, 0.2, 0))+
  geom_hline(yintercept = 0)+
  scale_fill_manual(values=c( blue1, red1), name=expression(paste(italic("in vivo")," response")))+
   ggtitle("Ibrutinib")+
  theme_CLL_ibr
# plot_ibr

remove(df_CLL_ibr, theme_CLL_ibr)

11.2 11.2 Venetoclax

# prepare CLL venetoclax dataset
df_CLL_ven<-data_all%>%mutate(Venetoclax=ifelse(grepl("Venetoclax", treatment_spec), 1, 0))
df_CLL_ven<-df_CLL_ven%>%filter(Drug=="Venetoclax")%>%filter(Venetoclax==1)%>%filter(diagnosis=="CLL")%>%filter(Response_not_evaluable==0)
df_CLL_ven$Resp_protocol<-factor(df_CLL_ven$Resp_protocol, levels=c("R"), 
                       labels=c("response        ")) # all patients had a R


plot_ven<-df_CLL_ven%>%ggplot(aes(x=reorder(patientID, -ES_AUC), y=-ES_AUC, fill=Resp_protocol)) +
  geom_bar(stat="identity", alpha=0.75)+
  scale_x_discrete(name="\nCLL patients")+
  scale_y_continuous(name=expression(paste("Effect size (",italic("ex vivo"),")")), limits=c(-0.75, 0), breaks=c(-0.75, -0.5, -0.25, 0), labels=c(0.75, 0.5, 0.25, 0))+
  geom_hline(yintercept = 0)+
  scale_fill_manual(values=c(blue1), name=expression(paste(italic("in vivo")," response")))+
  ggtitle("Venetoclax")+
  theme_figures_x_angle_45_wo_legend
 
# plot_ven

remove(df_CLL_ven)

11.3 11.3 Figure 5

design1<-"
AB
"

tp<-theme(plot.tag = element_text(size=22, face="bold"))

figure_CLL<-
  wrap_elements(plot_ven)+tp+
  wrap_elements(plot_ibr)+tp+
  plot_layout(design=design1)+
  plot_annotation(tag_levels ="A", title="Figure 5") & theme(title = element_text(size=18))

figure_CLL

# ggsave(path=(paste0(filepath,"04_Figures")), filename = "Figure_5_CLL_patients.pdf")
# remove(tp, plot_ibr, plot_ven, figure_CLL)

12 12 Quality assessment based on the type of sample material

plate_mean_normVal of inner DMSO controls by sample_type: BM, LN, PB

anova_res <- filter(screenData, Drug== "DMSO", !edge) %>% 
  group_by(plateID, material) %>% summarise(plate_mean_normVal = mean(normVal), .groups = "keep") %>% ungroup() %>%
  mutate(material = factor(material, levels = unique(material))) 

# anova_res %>% DT::datatable(width = 700)

Show ANOVA res

m <- lm(plate_mean_normVal ~ material, data = anova_res)
summary(m)
## 
## Call:
## lm(formula = plate_mean_normVal ~ material, data = anova_res)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.16263 -0.01092 -0.00358  0.00641  0.27585 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.00255    0.00313  320.77   <2e-16 ***
## materialBM  -0.00362    0.00763   -0.47     0.64    
## materialLN   0.00192    0.00925    0.21     0.84    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0348 on 162 degrees of freedom
## Multiple R-squared:  0.00186,    Adjusted R-squared:  -0.0105 
## F-statistic: 0.151 on 2 and 162 DF,  p-value: 0.86
remove(anova_res, m)

Show plate median sdev and range

filter(screenData, Drug== "DMSO", !edge) %>% group_by(plateID, patientID) %>% 
  summarise(plate_sdev = round(sd(normVal),4) , .groups = "keep") %>% ungroup() %>% 
  summarise(median_sdev = round(median(plate_sdev),4 ), 
            range_sdev = paste0( round(range(plate_sdev)[1],4) , 
                                 " - ", range(plate_sdev)[2] ), .groups = "keep") 
## # A tibble: 1 x 2
##   median_sdev range_sdev     
##         <dbl> <chr>          
## 1      0.0759 0.0293 - 0.6159

prepare data

plotTab <- select(screenData, plateID, plate_sdev, plate, material, frozen) %>% unique() %>%
  #change names of categories
  mutate(material = ifelse(material == "BM","Bone marrow",
                           ifelse(material == "LN","Lymph node", "Peripheral blood")), 
         frozen = ifelse(frozen == TRUE,"Frozen","Fresh"), 
         plate = ifelse(plate== "Acute","Myeloid","Lymphoid"))

show fig

p1<-ggplot(plotTab, aes(x = material, y = plate_sdev, fill = frozen, shape = plate)) +
  geom_point(aes(color= plate), alpha = 0, size=2.5) +
  geom_beeswarm(cex=1.5, size=2.5, alpha= 0.85, show.legend = F) +
  scale_shape_manual(values = c("Myeloid" = 25, "Lymphoid" = 21)) +
  #scale_color_manual(values = c("Fresh"="darkred","Frozen"="deepskyblue4")) +
  #scale_fill_manual(values = c("Fresh"="darkred","Frozen"="deepskyblue4")) +
  scale_color_manual(values = c("Fresh"=red2,"Frozen"=blue1)) +
  scale_fill_manual(values = c("Fresh"=red2,"Frozen"=blue1)) +
  geom_hline(yintercept = 0.3, linetype="22", col="gray30", size=0.7) +
  
  labs(color= "Material", shape = "Drug panel") +
    ylab("Standard deviation (SD)\n") +
    xlab("\nTumor sample origin") +
    guides(fill = "none", 
           color = guide_legend(override.aes = list(alpha = 1), reverse = T, nrow=2), 
           shape = guide_legend(override.aes = list(alpha = 1), nrow=2)) +
   theme_figures+
  theme(legend.position = "right")
p1  

# ggsave( path=(paste0(filepath,"04_Figures")), filename = "Figure_S1_Quality_assessment.pdf")


remove(plotTab)

12.1 12.1 Drug Drug Correlation

cor.mat <-
  select(data_all, patientID, Drug,  AUC) %>% 
  pivot_wider(names_from = Drug, values_from = AUC)%>% 
  column_to_rownames("patientID") %>% 
  rstatix::cor_mat()

cor.mat <- column_to_rownames(cor.mat, var="rowname")

breaks <- c(seq(-1, 1, length.out = 101)) %>% `names<-`(
     colorRampPalette(c("#003DA5", "#2055B0", "#406EBC", "#6086C7", "#809ED2", "#9FB6DD", "#BFCFE9", "#DFE7F4", "white", "#F4E0E7", "#E9C2CF", "#DEA3B6", "#D3849E", "#C76586", "#BC476E", "#B12855", "#A6093D"))(101))

pheatmap::pheatmap(cor.mat,
         breaks = breaks,
         color= names(breaks),
         border_color=NA,
         fontsize = 7,
         # filename=paste0(filepath,"04_Figures/Figure_S2.pdf"),
         width = 14,
         height = 14)

# ggsave( path=(paste0(filepath,"04_Figures")), filename = "Figure_S2.pdf")

# unique(data_all$patientID) %>% length()

13 13 Genotype - phenotype associations

13.1 13.1 FLT3-ITD ratio & FLT3 inhibitors

# prepare data
tmp<-data_all%>%
  # here all AML patients with an available information about FLT3-ITD mutation were considered, even if patients had no evaluable response they were included because here we have no       correlation with in-vivo response
  filter(diagnosis=="AML")%>%
  filter(!is.na(FLT3_ITD_ratio), FLT3_ITD_ratio!="na")%>%
  # only FLT3 inhibnitors are considered
  filter(Drug=="Crenolanib" | Drug=="Gilteritinib" | Drug=="Quizartinib"| Drug=="Sorafenib"| Drug=="Midostaurin hydrate")%>%
  mutate(Drug=as.character(Drug)) %>% 
  mutate(Drug=ifelse(Drug=="Midostaurin hydrate", "Midostaurin\nhydrate", Drug)) %>% 
  select(patientID, Target, diagnosis, FLT3_ITD_ratio, Drug_response_mean_low, Drug)

# correlation plot
plot_S2A<-tmp%>%ggplot(aes(x=Drug, y=Drug_response_mean_low, color=FLT3_ITD_ratio))+
  geom_beeswarm(size=3, alpha=1, cex=5) + 
  scale_color_gradientn(colors=c(blue1, red1), limits=c(0,1), name="FLT3-ITD ratio", breaks=c(0,0.5,1.0))+
  scale_y_continuous(name="Viability (AUC)")+
  theme_figures+
  stat_cor(aes(x=FLT3_ITD_ratio, y=Drug_response_mean_low), method="kendall",  cor.coef.name = c("tau"), label.sep="\n", label.x.npc=0.4, label.y.npc = 0.04, size=5)+
  facet_wrap(~ Drug, nrow=1, scales="free_x")+
  theme(axis.text.x = element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.x = element_blank())

plot_S2A

13.2 13.2 FLT3-TKD & FLT3 inhibitors

# prepare dataset
tmp<-data_all%>%
  # here all AML patients with an available information about FLT3-ITD mutation were considered, even if patients had no evaluable response they were included because here we have no       correlation with in-vivo response, only patients with FLT3-TKD mutation but not FLT3-ITD mutation are considered
  filter(diagnosis=="AML")%>%
  filter(!is.na(FLT3), FLT3!="na")%>%
  filter(FLT3_TKD==1 | FLT3==0)%>%
  # only FLT3 inhibnitors are considered
  filter(Drug=="Crenolanib" | Drug=="Gilteritinib" | Drug=="Quizartinib"| Drug=="Sorafenib"| Drug=="Midostaurin hydrate")%>%
  mutate(Drug=as.character(Drug)) %>% 
  mutate(Drug=ifelse(Drug=="Midostaurin hydrate", "Midostaurin\nhydrate", Drug)) %>% 
  select(patientID, Target, diagnosis, FLT3, Drug_response_mean_low, Drug)

# plot
plot_S2B<-ggplot(tmp, aes(x=Drug,y=Drug_response_mean_low,color=FLT3))+
  geom_beeswarm(size=2.5, alpha=1, cex = 1.8)+
  scale_color_manual(values=c(blue1, red1), name="", label=c("wildtype", "FLT3-TKD\nmutated"))+
  scale_x_discrete(name="")+
  scale_y_continuous(name="Viability (AUC)", limits=c(0, 1.4), breaks=c(0, 0.25, 0.5, 0.75, 1, 1.25))+
  geom_hline(aes(yintercept=1), colour="black", 
                linetype="dashed")  +
  theme_figures_x_angle_45+
  guides(color=guide_legend(title="FLT3", nrow=1))+
  theme(legend.position = "bottom", axis.title.x = element_blank() )

plot_S2B

remove(tmp)

13.3 13.3 IDH mutations and venetoclax

# prepare dataset
tmp<- data_all%>%
  # here all AML patients with an available information about IDH mutation were considered, even if patients had no evaluable response they were included because here we have no            correlation with in-vivo response
  filter(diagnosis=="AML")%>%
  filter(Drug=="Venetoclax")%>%
  select(patientID, IDH1, IDH2, Drug, Pathway, AUC)%>%
  mutate(IDH=case_when(IDH1==1 ~ "IDH1\nmutated", IDH2==1 ~ "IDH2\nmutated", IDH1==0 & IDH2==0 ~ "wildtype", IDH1=="na" & IDH2=="na"~ "NA" ))%>%
  filter(IDH!="NA")

# plot
plot_S2C<-ggplot(tmp, aes(x=IDH,y=AUC,color=IDH))+
# geom_boxplot()+
  geom_boxplot(data=summarise(group_by(tmp, IDH), mean=mean(AUC), sd=sd(AUC), .groups = "keep"),
               aes(x=IDH,color=IDH,lower=mean-sd,upper=mean+sd,middle=mean,ymin=mean-sd,ymax=mean+sd),
               stat="identity", inherit.aes=FALSE)+
  geom_beeswarm(size=3, cex=3)+
 scale_color_manual(values=c(red1, gray1, blue1))+
  scale_y_continuous(name="Viability (AUC)", limits=c(0, 1.2), breaks=c(0, 0.25, 0.5, 0.75, 1))+
  geom_hline(aes(yintercept=1), colour="black", 
                linetype="dashed")  +
  ggpubr::stat_compare_means(method = "t.test",paired=F,  size=starsize, comparisons = list( c(1,3)), method.args = list(       alternative = c("two.sided")))+
  theme_figures_wo_legend+
  xlab("")

plot_S2C

remove(tmp)

13.4 13.4 Nutlin-3a and TP53

# prepare dataset
tmp<-data_all%>%
  filter(!is.na(TP53), TP53!="na")%>%group_by(patientID,  TP53)%>%
  filter(Drug=="Nutlin-3a")%>%mutate(TP53=case_when(TP53==1 ~ "TP53\nmutated", TP53==0 ~ "wildtype"))

# plot
plot_S2D<-ggplot(tmp, aes(x= TP53,
             y=AUC, 
             color=TP53))+
  geom_boxplot(data=summarise(group_by(tmp, TP53), mean=mean(AUC), sd=sd(AUC), .groups = "keep"),
               aes(x=TP53,color=TP53,lower=mean-sd,upper=mean+sd,middle=mean,ymin=mean-sd,ymax=mean+sd),
               stat="identity", inherit.aes=FALSE)+
  geom_beeswarm(size=3, cex=3)+
  scale_y_continuous(name="Viability (AUC)", limits=c(0, 1.2), breaks=c(0, 0.25, 0.5, 0.75, 1))+
  geom_hline(aes(yintercept=1), colour="black", 
                linetype="dashed")  +
  scale_x_discrete(name="")+
  scale_color_manual(values=c(red1, blue1))+
  ggpubr::stat_compare_means(method = "t.test",paired=F,  size=starsize, comparisons = list( c(1,2)), method.args = list(       alternative = c("two.sided")))+
  theme_figures_wo_legend
plot_S2D

remove(tmp)

13.5 13.5 Plot Figure S2

design1<-"
AA
BC
BD
"
tp<-theme(plot.tag = element_text(size=22, face="bold"))


 wrap_elements(plot_S2A)+tp+
 wrap_elements( plot_S2B)+tp+
 wrap_elements( plot_S2C)+tp+
 wrap_elements( plot_S2D)+tp+
 plot_layout(design=design1, heights = c(5,3, 3))+
  plot_annotation(tag_levels = "A")

# ggsave( path=(paste0(filepath,"04_Figures")), filename = "Figure_S3.pdf")
  
#ggsave(width = 16, height =24, path=(paste0(filepath,"04_Figures")), filename = "Figure_S2.pdf")
# remove(tp, plot_S2A, plot_ S2B, plot_S2C, plot_S2D, figure_S2, plot_blank)

14 14 Figure S3: Boxplots for all significant drugs: R vs SD vs PD

prepare data

#list of drugs that passed 0.05 (R vs PD). Arrange by pathway alphabetical (as fig 3B) and 
#then by smallest pval.
example_drugs <- filter(df_p_drugAUC, pval <= 0.05) %>% arrange( Pathway , pval ) %>% 
  select(name) %>% unique()

#get the df
plotTab <- filter(df_chemo_bckup, name %in% example_drugs$name) %>% 
  select(-concentration, -concIndex, -normVal) %>% unique() %>%
  mutate(plotTxt = paste0(Pathway," - ", name), 
         box_order = ifelse(response == "PD",0,ifelse(response == "SD", 1,2)) ) %>% 
  arrange( box_order ) %>% mutate(response = factor(response, levels=unique(response)))

remove(df_p_drugAUC)

boxplots

#make plot list
plotList <- lapply( example_drugs$name ,function(drugname){
  
  p <- ggplot(plotTab[plotTab$name==drugname,], aes(x=response, y=AUC)) + 
    geom_boxplot(aes(fill = response), width = 0.8, size=0.3, alpha = 0.85, 
                 color = "black", outlier.shape = NA) +
    scale_fill_manual(values = c("PD"=red1,"SD"=gray1, 
                                 "R"=blue1) ) +
    geom_beeswarm(aes(color = response),cex=2, dodge.width=0.8, shape=21, 
                  fill="slategray4", size=2, alpha= 0.7) +
    scale_color_manual(values = c("black", "black", "black") ) + 
    geom_hline(yintercept = 1, linetype="22", col="gray30", size=0.7) +
    
    coord_cartesian(ylim = c(0,1.5)) +
    
    labs(title= drugname ) +
    xlab(expression(paste(italic("in vivo")," response"))) +
    ylab("Viability (AUC)") +
    guides(color="none", fill="none") +
    theme_classic()+
    theme( axis.line = element_blank(), 
           panel.border = element_rect(size=1, fill = NA),
            plot.title = element_text(face='bold', hjust = 0.5,  size = 14))
  return( p )
})

See examples

cowplot::plot_grid(plotlist = plotList[1:8], ncol=2)

save in pdf

makepdf(plotList,  "../../inst/Figure_S4_boxplots.pdf", ncol =2,nrow = 4,
                 figNum = 8, width =8, height=12)

remove(plotList, plotTab, example_drugs, df_chemo_bckup)

15 15 Effect of Tumor Infiltration

Prepare data

df_tinf <- filter(screenData, plate_sdev < 0.3, Drug != "DMSO", 
                  Pathway == "Chemotherapy", chemo_pat == T) %>% 
  select(patientID,  diagnosis, treatment_type, chemo_pat, resp_protocol,
         t_infiltration, plateID, Drug, concentration, normVal) %>% 
  #one viability per drug conc per patient 
  #(average, for drugs in both plates and replicate plates)
  group_by(patientID, Drug, concentration) %>% mutate(normVal = mean(normVal)) %>% 
  select(-plateID) %>% unique() %>% 
  #AUC drug pat
  group_by(patientID, Drug) %>% mutate(AUC = calcAUC(normVal, concentration)) %>% 
  select(-concentration, -normVal) %>% unique() %>%
  #one viability per pat: mean AUC pat
  group_by(patientID) %>% mutate(AUC_pat = mean(AUC)) %>% select(-Drug, -AUC) %>% unique() %>%
  ungroup() %>%
  mutate(t_infiltration = as.numeric(t_infiltration)) %>% #group_by(resp_protocol) %>%
  #correlation test
  mutate(corr_resp = round(cor(t_infiltration, AUC_pat, method = c("pearson")),2), 
         corr_pval = round(cor.test(x = t_infiltration, y = AUC_pat, method = c("pearson"))$p.value,2), 
         txt = paste0("R = ",corr_resp, "\np = ", corr_pval)) 

Plot New Plot not subsetting by response group

p1<-ggplot(df_tinf, aes(x=t_infiltration, y=AUC_pat)) + 
  geom_point(size=2.5, show.legend = T, alpha=0.85) +
  
  geom_smooth(method = "lm", se=F, show.legend = F, formula='y ~ x')+
  
  annotate("text", x = 60, y = 0.9,
  label = c(unique(df_tinf[,]$txt)) ,
  color=c("black"), size=5) +
  
  coord_cartesian(xlim = c(50,100), clip = "off") +
  xlab("\nTumor cell infiltration (in %)") + ylab("Viability (AUC)\n") +
  theme_figures+
  theme(aspect.ratio = 1)

# ggsave( path=(paste0(filepath,"04_Figures")), filename = "Figure_S4_Tumor_cell_infiltration.pdf")

p1

plotTab<-filter(screenData, plate_sdev < 0.3, Drug != "DMSO", 
                  Pathway == "Chemotherapy", chemo_pat == T) %>% 
  select(patientID,  diagnosis, treatment_type, chemo_pat, resp_protocol,
         t_infiltration) %>% 
  group_by(patientID,  diagnosis, treatment_type, chemo_pat, resp_protocol,
         t_infiltration) %>% 
  summarize(.groups = "keep") %>% 
  filter(resp_protocol%in%c("R","PD")) 
  

plotTabmeanboxplot<-plotTab %>% 
  group_by(resp_protocol) %>% 
  summarise(mean=mean(t_infiltration), sd=sd(t_infiltration))

p2<-  ggplot(plotTab, aes(x=resp_protocol,color=resp_protocol,y=t_infiltration))+
  geom_boxplot(data=plotTabmeanboxplot, aes(x=resp_protocol,color=resp_protocol,lower=mean-sd,upper=mean+sd,middle=mean,ymin=mean-sd,ymax=mean+sd),stat="identity", inherit.aes 
=FALSE)+
  # geom_boxplot()+
  geom_beeswarm(aes())+
  stat_compare_means(comparisons = list(c(1,2)), method.args = list(       alternative = c("two.sided")))+
  theme_figures+
  guides(color="none")+
  ylab("\nTumor cell infiltration (in %)") +
    xlab(expression(paste(italic("in vivo")," response"))) +
    scale_color_manual(values = c("R" = blue1, "PD"=red2))


p2

#eliminate conc-specific rows/cols          
df_p_drugAUC <- select(df_chemo, 
                       -concentration, -concIndex, -normVal, -AUC_pathway_patient) %>% unique() %>%
  #add tumor infiltration variable
  mutate(t_infiltration = as.numeric(screenData[match(patientID, screenData$patientID),]$t_infiltration) ) %>%
  #"AUC" to "metric" for functions
  mutate(metric = AUC) %>% 
  #perform tests 
  group_by(name) %>% nest() %>%
         #t-test
  mutate(pval_t = map(data, t_test), 
         #multivariate regression
         pval_mreg = map(data, m_reg)) %>% 
  unnest( c(data, pval_t, pval_mreg) ) %>% ungroup() %>%
  #eliminate patient data
  select(Pathway, name, pval_t, pval_mreg) %>% unique() %>%
  #add annotations
  mutate(label_col = ifelse(pval_t > 0.05 & pval_mreg > 0.05, "01both_not_sig", 
                              ifelse( pval_t < 0.05 & pval_mreg < 0.05, "02both_sig", 
                                      ifelse(pval_t < 0.05 & pval_mreg > 0.05, "03sig_t_only", 
                                             ifelse( pval_t > 0.05 & pval_mreg < 0.05, "04sig_mreg_only","") ) ) )) %>%
  #plotting order of points
  arrange(label_col)
  
# df_p_drugAUC %>%
#   mutate_if(is.numeric, formatC, digits=2) %>% DT::datatable(width = 700)
p3<-ggplot(df_p_drugAUC, aes(x=-log10(pval_t) , y=-log10(pval_mreg), fill = label_col, color = label_col )) + 
  geom_point(size=3, shape = 21, alpha = 0.8, show.legend = T) + 
  
  scale_fill_manual(values = c("01both_not_sig"="#E18727FF","02both_sig"="deepskyblue4", "03sig_t_only"="darkred"), 
                     labels = c("01both_not_sig"="p>0.05 in both models", 
                                "02both_sig"="p<0.05 in both models",
                                "03sig_t_only"="p<0.05 without accounting\nfor tumor infiltration")) +
  scale_color_manual(values = c("01both_not_sig"="#E18727FF","02both_sig"="deepskyblue4", "03sig_t_only"="darkred"), 
                     labels = c("01both_not_sig"="p>0.05 in both models", 
                                "02both_sig"="p<0.05 in both models",
                                "03sig_t_only"="p<0.05 without accounting\nfor tumor infiltration")) +
  
  geom_abline(intercept = 0, slope = 1, linetype="22", col="gray30", size=0.7) +

  xlab(expression(-log[10]*'('*italic(P)*'), tumor infiltration not considered')) +
  ylab(expression(-log[10]*'('*italic(P)*'), accounting for tumor infiltration')) +
  labs(fill = "Statistical significance",
       color = "Statistical significance") +
  theme_figures+
  theme(legend.position = "right", aspect.ratio = 1)+
  guides(fill=guide_legend(nrow=3))

p3

tp<-theme(plot.tag = element_text(size=22, face="bold"))

design2<-"
AB
CC"

wrap_elements(p1)+tp+
  wrap_elements(p2)+tp+
  wrap_elements(p3)+tp+
  plot_layout(design=design2, heights = c(1,1), widths = c(2,1))+
  plot_annotation(tag_levels = "A")
## Warning in wilcox.test.default(c(88, 70, 59, 80, 88, 90, 74), c(83, 93, : cannot
## compute exact p-value with ties

# ggsave( path=(paste0(filepath,"04_Figures")), filename = "Figure_S5.pdf")

16 16 Session info

sessionInfo()
## R version 4.1.3 (2022-03-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
## [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
## [5] LC_TIME=German_Germany.1252    
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] Cairo_1.5-15       writexl_1.4.0      corrplot_0.92      glmnet_4.1-3      
##  [5] Matrix_1.4-1       table1_1.4.2       patchwork_1.1.1    forestplot_2.0.1  
##  [9] checkmate_2.0.0    magrittr_2.0.3     forestmodel_0.6.2  maxstat_0.7-25    
## [13] survival_3.3-1     survminer_0.4.9    ggpubr_0.4.0       dplyr_1.0.8       
## [17] forcats_0.5.1      stringr_1.4.0      purrr_0.3.4        readr_2.1.2       
## [21] tidyr_1.2.0        tibble_3.1.6       tidyverse_1.3.1    drc_3.0-1         
## [25] MASS_7.3-56        gridExtra_2.3      cowplot_1.1.1      ggbeeswarm_0.6.0  
## [29] gghighlight_0.3.2  ggrepel_0.9.1      RColorBrewer_1.1-3 ggpmisc_0.4.6     
## [33] ggpp_0.4.4         ggplot2_3.3.5      knitr_1.38         BiocStyle_2.22.0  
## 
## loaded via a namespace (and not attached):
##  [1] TH.data_1.1-1         colorspace_2.0-3      ggsignif_0.6.3       
##  [4] ellipsis_0.3.2        markdown_1.1          fs_1.5.2             
##  [7] gridtext_0.1.4        ggtext_0.1.1          rstudioapi_0.13      
## [10] farver_2.1.0          MatrixModels_0.5-0    DT_0.22              
## [13] fansi_1.0.3           mvtnorm_1.1-3         lubridate_1.8.0      
## [16] xml2_1.3.3            codetools_0.2-18      splines_4.1.3        
## [19] Formula_1.2-4         jsonlite_1.8.0        broom_0.7.12         
## [22] km.ci_0.5-6           cluster_2.1.3         dbplyr_2.1.1         
## [25] pheatmap_1.0.12       BiocManager_1.30.16   compiler_4.1.3       
## [28] httr_1.4.2            backports_1.4.1       assertthat_0.2.1     
## [31] fastmap_1.1.0         cli_3.2.0             htmltools_0.5.2      
## [34] quantreg_5.88         tools_4.1.3           gtable_0.3.0         
## [37] glue_1.6.2            Rcpp_1.0.8.3          carData_3.0-5        
## [40] cellranger_1.1.0      jquerylib_0.1.4       vctrs_0.4.0          
## [43] nlme_3.1-157          crosstalk_1.2.0       iterators_1.0.14     
## [46] exactRankTests_0.8-34 xfun_0.30             rvest_1.0.2          
## [49] lifecycle_1.0.1       gtools_3.9.2          rstatix_0.7.0        
## [52] zoo_1.8-9             scales_1.1.1          hms_1.1.1            
## [55] sandwich_3.0-1        SparseM_1.81          yaml_2.3.5           
## [58] KMsurv_0.1-5          sass_0.4.1            stringi_1.7.6        
## [61] highr_0.9             foreach_1.5.2         plotrix_3.8-2        
## [64] shape_1.4.6           rlang_1.0.2           pkgconfig_2.0.3      
## [67] evaluate_0.15         lattice_0.20-45       htmlwidgets_1.5.4    
## [70] labeling_0.4.2        tidyselect_1.1.2      bookdown_0.25        
## [73] R6_2.5.1              magick_2.7.3          generics_0.1.2       
## [76] multcomp_1.4-19       DBI_1.1.2             mgcv_1.8-40          
## [79] pillar_1.7.0          haven_2.4.3           withr_2.5.0          
## [82] abind_1.4-5           pamr_1.56.1           modelr_0.1.8         
## [85] crayon_1.5.1          car_3.0-12            survMisc_0.5.6       
## [88] utf8_1.2.2            tzdb_0.3.0            rmarkdown_2.13       
## [91] readxl_1.4.0          data.table_1.14.2     reprex_2.0.1         
## [94] digest_0.6.29         xtable_1.8-4          munsell_0.5.0        
## [97] beeswarm_0.4.0        vipor_0.4.5           bslib_0.3.1